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Finite  Difference  Modelling  of 
Scattering  by  Objects  in  the  Seabed 

John  A.  Fawcett  and  Joris  L.T.  Grimbergen 


Executive  Summary; 

A  better  understanding  of  the  scattering  of  acoustic  energy  from  elastic  objects 
on  or  under  the  seabed  will  lead  to  a  significant  improvement  in  our  ability  to 
detect  and  classify  minelike  objects.  Modelling  of  such  scattering  should  include 
the  effects  of  the  seabed  itself  since  mines  may  be  buried  to  some  degree.  The 
Finite  Difference  method,  described  in  this  report,  allows  acoustic  scattering 
from  mines  in  complex  bathymetric  conditions  to  be  modelled.  This  gives 
insight  into  the  physical  mechanisms  and  environmental  parameters  which  are 
important  to  the  scattering  of  energy  from  mines  and  will  act  as  a  benchmark 
for  faster,  but  more  approximate,  models. 

This  report  describes  some  of  the  theory  and  implementation  issues  concerned 
with  Finite  Difference  modelling.  For  an  illustration  of  the  value  of  this  method, 
the  results  of  computations  of  scattering  from  buried  and  partially  buried  cylin¬ 
ders  for  smooth  and  rough  seabed  conditions  are  presented.  Numerical  exam¬ 
ples  are  also  given,  illustrating  the  accuracy  of  the  method  for  problems,  such 
as  scattering  from  aluminum  cylinders,  where  analytic  solutions  are  known.  It 
is  clear  many  questions  regarding  the  effects  of  burial  and  bathymetry  on  mine 
scattering  are  answered  by  the  FD  model. 
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Finite  Difference  Modelling  of 
Scattering  by  Objects  in  the  Seabed 

John  A.  Fawcett  and  Joris  L.T.  Grimbergen 


Abstract:  In  this  report  we  describe  the  theory  and  some  implementation 

issues  of  a  finite  difference  code  used  at  SACLANT  Centre.  In  particular,  we 
consider  the  modelling  of  attenuation  and  the  excitation  of  a  remote  incident 
field  by  using  Huy  gen’s  sources.  A  series  of  comparisons  of  finite  difference 
results  with  analytical  results  is  performed.  The  report  concludes  with  a  series 
of  computations  of  scattering  of  a  generalized  plane  wave  from  a  buried  cylinder 
where  the  transmitted  portion  of  the  generalized  plane  wave  is  evanescent.  An 
example  of  time-reversed  propagation  of  a  scattered  field  is  also  given. 

Keywords:  Finite  difference  modelling,  attenuation,  Huygen’s  sources 
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1 

Introduction 


Finite  Difference  (FD)  modelling  can  compute  the  time-domain  solution  to  complex 
elastic  scattering  problems  (for  example,  [1,  2,  3,  4]  ).  The  FD  solutions  include 
all  angles  of  energy  propagation  -  both  forward  and  backward  propagating  energy. 
The  method  is  very  general  in  terms  of  the  environment  and  the  scattering  objects 
and  surfaces.  There  is  no  need  to  make  assumptions  about  range-independence,  the 
shapes  of  objects,  etc.  The  chief  disadvantage  of  the  method  is  that  it  becomes 
computationally  intensive  for  grid  sizes  which  are  more  than  100  wavelengths  in  a 
dimension  (for  three-dimensional  modelling  this  problem  is  even  more  acute).  Thus 
for  problems  where  the  source  of  incident  energy  is  distant  from  the  scattering  region 
it  is  not  possible  to  model  the  source  region.  In  order  to  resolve  this  shortcoming,  we 
have  implemented  a  technique  for  bringing  remote  incident  fields  into  the  numerical 
grid.  This  means  that  in  the  case  where  a  large  portion  of  the  waveguide  is  simple  in 
structure,  we  can  model  the  propagation  of  the  incident  field  up  to  the  boundary  of 
the  grid  with  a  more  efficient  technique.  Similarly,  once  the  scattered  field  has  been 
computed  by  the  FD  method,  it  should  be  possible  to  extrapolate  the  scattered  field 
to  remote  receivers. 

In  the  Spring  of  1995,  SACLANTCEN  obtained  the  visco-elastic  finite  difference 
code  developed  by  J.  Robertsson,  J.  Blanch,  and  W.  Symes  [5,  6].  This  method  in¬ 
cludes  additional  variables  and  parameters  which  allow  for  the  modelling  of  spatially- 
varying  compressional  and  shear  attenuations.  We  shall  refer  to  this  code  as  the 
RBS  code  for  the  remainder  of  this  report.  We  implemented  the  basic  code  into  a 
MATLAB  [7]  package  for  FD  modeUing.  This  package  provides  subroutines  for  the 
construction  of  the  fields  required  by  the  FD  code.  This  includes  the  construction 
of  stress  and  strain  relaxation  time  grids  for  desired  Q  (the  number  of  wavelengths 
over  which  the  amplitude  decays  by  e“^)  grids  for  the  compressional  and  shear 
waves.  It  also  includes  a  subroutine  for  the  analytic  computation  of  the  generalized 
plane- wave  incident  field  along  a  boundary  of  the  finite  difference  grid  for  a  2  half¬ 
space  environment.  The  incident  field,  in  this  case,  consists  of  a  coherent  sum  of 
direct,  refiected,  and  transmitted  components.  This  field  is  used  to  excite  sources 
along  this  boundary  in  the  FD  computations,  effectively  bringing  the  incident  field 
into  the  grid.  The  FD  code  is  fourth-order  accurate  spatially  and  uses  a  staggered 
grid  formulation  [4].  It  is  second  order  in  time.  In  many  instances,  20  grid  points 
per  wavelength  give  very  good  results  with  this  code  (  for  a  broadband  source,  this 
corresponds  to  fewer  points  per  wavelength  for  frequencies  higher  than  the  centre 
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value). 

Below,  we  will  describe  in  more  detail  some  of  the  features  of  the  FD  code;  in 
particular,  we  will  describe  some  of  the  details  of  the  attenuation  modelling,  the 
boundary  absorbing  layer,  and  the  incident  field  modelling.  We  will  then  present  the 
results  from  a  sequence  of  benchmark  computations  which  test  the  attenuation  and 
the  elastic  scattering  capabilities  of  the  code.  In  order  to  reduce  computation  time, 
the  numerical  examples  of  this  report  were  run  using  a  FORTRAN  implementation 
of  the  code  rather  than  the  MATLAB  package. 
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2 

Theory 


In  this  section  we  describe  some  of  the  theoretical  and  implementation  issues  con¬ 
nected  with  the  visco-elastic  finite  difference  modelling. 


2.1  The  basic  mode! 


The  basic  set  of  visco-elastic  differential  equations  solved  by  the  RBS  code  is 
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where  denotes  the  components  of  the  symmetric  stress  tensor,  Vi  is  the  velocity 
vector,  p  is  the  shear  relaxation  modulus,  and  7/  is  the  compressional  relaxation 
modulus,  corresponding  to  A  +  2/x  in  the  standard  elastic  problem.  The  variables 
Vij  are  memory  variables  [8]  which  are  introduced  for  the  visco-elastic  modelling  in 
order  to  avoid  the  explicit  computation  of  time-convolution  terms.  For  the  simple 
attenuation  model  considered  in  this  report,  three  memory  variables  are  required. 
Associated  with  the  attenuation  modelling  are  the  parameters  the  compressional 
strain  relaxation  time,  the  shear  strain  relaxation  time,  and  T(j  the  stress  relax¬ 
ation  time  which  for  our  model  is  the  same  for  compressional  and  shear  waves.  In 
the  limit  of  no  attenuation,  and  become  unity  and  Vij  becomes  zero. 
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The  above  equations  then  reduce  to  the  standard  equations  of  elasticity.  In  the 
following  subsection,  we  discuss  in  more  detail  the  modelling  of  attenuation. 

As  with  all  FD  codes  the  differential  equations  are  discretized  in  both  space  and  time 
to  yield  a  set  of  discrete  equations.  The  RBS  code  uses  a  discretization  of  Eq.(l) 
which  is  fourth  order  accurate  in  space  and  second  order  in  time.  In  order  to  obtain 
accurate  answers  using  the  FD  code  we  find  in  practice  that  we  require  10-20  spatial 
grid  points  for  the  dominant  wavelength  of  the  problem.  This  criterion  is  dependent 
upon  the  frequency  bandwidth  of  the  source;  if,  for  example,  there  is  significant 
energy  at  frequencies  twice  the  main  frequency,  these  wavelengths  must  also  be 
accurately  modelled.  The  required  spatial  discretization  of  a  wavelength  indicates 
the  problem  which  arises  when  dealing  with  low-velocity  zones;  the  wavelengths  are 
small  in  these  regions  and  hence  these  zones  require  a  small  spatial  step.  The  code 
of  this  report  uses  uniform  spacing,  so  that,  in  fact,  a  small  spacing  must  be  used 
for  the  entire  grid  in  this  case.  These  low  velocity  cases  may  arise  when  we  wish  to 
model  an  object  fiUed  with  air  (cp  ^  340m/ s)  or  when  we  wish  to  accurately  model 
the  effects  of  shear.  Shear  velocities  may  often  be  only  a  few  hundred  m/s. 

The  RBS  is  an  explicit  FD  code  and  hence  the  time  step  At  must  satisfy  a  constraint 
of  the  form  At  <  aAx/cjy^ax  where  Ax  is  the  spatial  step  size  if  the  solution  is  to 
be  stable  with  respect  to  time.  This  has  two  main  implications  for  the  modelling; 
(1)  if  one  decreases  the  spatial  grid  size,  then  it  is  necessary  to  decrease  the  time 
step  accordingly  (2)  a  zone  of  high  velocity  will  require  small  time  steps  and  since 
the  time  step  is  the  same  for  the  entire  grid,  this  may  force  the  time  step  to  be 
inordinately  small  in  other  regions  of  the  grid.  For  the  RBS  code  the  maximum 
value  of  a  is  approximately  0.606  [9]. 


2.2  Attenuation  Modelling 

We  will  give  a  brief  outline  of  the  theory  of  visco-elastic  modeUing  to  give  an  under¬ 
standing  of  some  of  the  parameters  and  implementation  issues  involved.  Here,  we 
follow  the  notation  of  the  appendix  in  [6].  The  constitutive  relation  for  a  visco-elastic 
medium  can  be  modelled  as  being  time-dependent, 

O^ij  —  A  ^  ^kk^ij  T  2il^  ^  €ij  (^) 

where  aij  is  the  stress  tensor  of  Eqs.(l),  eij  is  the  strain  tensor,  A  and  M  correspond 
to  the  Lame  constants  A  and  /x,  and  Sij  denotes  the  Kronecker  delta.  Equation  (2) 
can  be  rewritten  in  terms  of  compressional  F  and  shear  parameters  M  as 

aij  =  (f  —  2M)  *  Ckk^ij  -f-  2M  *  €ij.  (3) 

The  form  of  Eqs.(2)  and  (3)  is  similar  to  the  standard  constitutive  relations  except 
that  now  the  parameters  A,  F,  and  M  are  time-dependent  and  multiplications  have 
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become  time- convolutions.  A  useful  model  for  the  time  history  of  F  and  M  is 

L  /  _p 


r(/)  =  Tr 


U{t) 


(4) 


and 


Mit)  =  Mr 


L 


^=1 


Tal) 


V(t) 


(5) 


where  U{t)  is  Heaviside’s  function,  and  are  the  compressional  and  shear  strain 
relaxation  times  for  the  £’th  term  in  Eqs.(4'5),  and  are  the  corresponding  stress 
relaxation  times.  Following  [6]  we  have  chosen  the  parameters  to  be  the  same  for 
both  the  compressional  and  shear  functions  and  hence  we  do  not  use  a  superscript 
for  this  parameter.  In  the  frequency  domain  the  complex  compressional  and  shear 
moduli  are  given  by 


11 

'immt)) 

d 

M{u)  =  F 

(6) 

(7) 


where  F  denotes  the  Fourier  Transform.  In  order  to  simplify  our  analysis  we  will 
consider  -£  =  1  in  Eqs.  (4-5)  and  wiU  no  longer  use  the  subscript  i  in  our  notation.  We 
will  consider  only  the  compressional  term  r(a;).  The  analysis  of  the  shear  function 
proceeds  along  identical  lines.  Using  Eq.(4)  in  Eq.(6)  yields 


f(cu) 


1  +  ioJT^ 

1  -1-  iuT^  ' 


(8) 


In  terms  of  r(u;)  the  quality  factor  Q  can  be  written 


n(,  A  = 


(9) 


where  and  S  denote  the  real  and  imaginary  parts  of  a  complex  number  respectively. 
This  factor  is  the  number  of  wavelengths  a  harmonic  plane-wave  must  propagate 
before  its  amplitude  has  been  decreased  by  a  factor  of  exp(— tt).  The  phase  velocity 
c(u;)  for  each  Fourier  component  of  the  wavefield  can  also  be  computed  from  T{u) 
and  is  given  by 


Te  Mr 
ra  P 


1  + 


(10) 


From  Eq.(lO)  it  can  be  seen  that  there  is  frequency  dispersion  in  this  visco-elastic 
medium.  In  order  to  have  a  causal  system  it  is,  in  fact,  necessary  to  have  some 
frequency  dispersion[10]. 


Introducing 


(11) 
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into  Eq.(9),  we  can  write 


l+a;2r2(l  +  r^)' 


For  reasonable  values  of  Q,  < <  1  and  hence 


Q~\^) 


UJTfjT^ 

1  +  ' 


(12) 


(13) 


In  order  to  determine  optimal  values  of  To-  and  to  approximate  a  constant  Q  over  a 
given  frequency  band,  [6]  gives  a  simple,  but  effective  algorithm  based  upon  Eq.(13). 
Basically  controls  the  amplitude  of  and  controls  the  frequency  offset  of 

the  curve.  The  parameter  r^r  is  set  to  1/(2%  fc)  where  fc  is  the  centre  frequency  of  the 
frequency  interval  of  interest  and  is  determined  by  a  simple  linear  optimization 
algorithm. 


The  same  type  of  algorithm  can  also  be  used  in  the  case  of  more  than  one  (i.e., 
^  >  1  in  Eq.(4))  attenuation  mechanism:  are  set  for  specified  frequencies  over  the 

bandwidth  and  a  value  of  is  determined  to  yield  a  good  approximation  to  constant 
Q  over  the  entire  frequency  band.  By  taking  several  attenuation  mechanisms,  it  is 
possible  to  obtain  an  excellent  approximation  to  a  constant  Q.  In  the  numerical  code 
of  this  report  only  one  mechanism  is  used.  Although  this  does  not  model  a  constant 
Q  over  the  entire  band  of  interest,  we  will  see  in  the  numerical  examples  that,  in 
practice,  it  does  a  good  job  in  modelling  the  attenuation  and  dispersion  effects  of 
visco-elasticity. 

In  Fig.  1  we  plot  Eq.(9)  as  a  function  of  frequency.  We  have  determined  the  optimal 
values  of  for  the  frequency  interval  [100,3000]  Hz  for  constant  Q  values  of  20 
(which  corresponds  to  1.36  dB/A)  and  50  (which  corresponds  to  0.54  dB/A).  The  Q 
for  the  single  mechanism  model  is  approximately  constant  in  the  interval  [750,3000] 
Hz.  For  frequencies  higher  or  lower  than  these  values,  the  Q  values  produced  by 
the  model  are  too  high.  If  the  frequency  content  of  the  source  is  concentrated  in 
the  interval  where  the  constant  Q  approximation  is  good,  then  we  expect  the  single¬ 
mechanism  Q  model  to  be  effective. 

In  Fig.  2  we  show  the  velocity /frequency  curve  computed  from  Eq.(lO)  for  the  com¬ 
puted  values  of  r^r  and  for  the  Q  =  20  (solid)  and  Q  =  50  (dashed)  cases.  The 
phase  velocity  varies  from  about  1500  m/s  at  /=:300  Hz  to  1575  m/s  at  /=5000  Hz 
for  the  Q  =  20  case  and  the  variation  is  almost  linear  between  500  and  1500  Hz. 
The  velocity  variation  is  much  less  for  the  Q  =  50  case;  there  is  only  a  variation  of 
30  m/s  over  the  entire  frequency  range. 

Based  on  the  theoretical  analysis  of  Futterman  [11]  and  the  experimental  work  of 
Wuenschel  [12],  the  following  Q/dispersion  model  for  real  visco-elastic  materials  has 
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Figure  1:  Computed  Q  as  a  function  of  frequency,  approximating  Q=20  (solid  line), 
Q=50  (dashed  line) 


Figure  2:  Computed  phase  velocity  as  a  function  of  frequency  for  modelled  Q=20 
(solid  line),  Q=50  (dashed  line) 
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been  suggested 
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27rQo 
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1-1 
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(14) 

(15) 


There  are  three  parameters  in  Eqs.  (14)  and  (15)  to  choose,  namely  u;o,  cq  and 


Figure  3:  Computed  phase  velocity  as  a  function  of  frequency  for  Q=20  using  Eq.(lO) 
(solid  line)  and  Eq,(14)  (dashed  line) 

Qo-  The  parameter  Uo  is  a  frequency  taken  to  be  weU  below  the  frequency  band 
of  interest,  cq  is  a  velocity  chosen  to  yield  a  desired  phase  velocity  at  a  specific 
frequency  and  Qo  is  some  constant  Q- value.  In  Figs.  3  and  4  we  show  the  curves 
computed  from  Eqs.  (14)  and  (15)  with  the  corresponding  curves  for  Q{oj)  and  c(u?) 
in  Figs.l  and  2  for  a  value  of  Q  ~  20.  We  have  used  Qo  =  20,  uq  ==  OAHz,  and  co 
is  chosen  such  that  c{2Q0Hz)  =  lA85m/s  in  Eq.  (14)  and  (15). 

From  Fig.  3  we  can  see  that  with  the  appropriate  choice  of  parameters  the  curve  com¬ 
puted  from  Eq.  (14)  agrees  well  with  the  curve  computed  for  the  single-mechanism 
model  except  at  the  low  frequencies.  The  value  of  Q{uj)  computed  from  Eq.  (15) 
is  essentially  constant,  with  a  value  slightly  less  than  20,  over  the  frequency  range 
shown. 

We  have  analyzed  the  characteristics  of  the  single-mechanism  model  in  the  frequency 
domain.  However,  the  Finite  Difference  code  is  implemented  in  the  time  domain. 
From  Eq.  (3)  this  would  seem  to  require  the  computation  of  time-convolution  terms 
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Figure  4:  Computed  Q  as  a  function  of  frequency  for  Q=20  using  Eq,(9)  (solid  line) 
and  Eq.(15)  (dashed  line) 


at  each  spatial  grid  point.  Fortunately,  due  to  the  form  o{T{t)  and  M{t)  in  Eqs.  (4) 
and  (5)  each  of  the  convolution  expressions  for  a^z  and  a^x  can  be  expressed 
in  terms  of  a  time-dependent  term  and  time-dependent  variables  r^x^,  t'zz  and  r^z 
respectively  which  satisfy  straightforward  differential  equations.  Thus,  instead  of  ex¬ 
plicitly  computing  time  convolutions,  the  differential  equations  for  the  convolutions 
are  updated  at  the  same  time  steps  as  the  standard,  elastic  differential  equations. 
For  the  single-mechanism  model,  three  extra  differential  equations  are  required  (at 
each  spatial  point). 

In  the  MATLAB  implementation  of  the  code,  a  module  has  been  written  which 
takes  user-specified  grids  of  Q-values  and  sound  speeds  and  converts  this  into  grids 
of  Te  and  Tfj  values.  These  grids  are  required  by  the  FD  code  for  the  attenuation 
modelling. 


2.3  Boundary  attenuation 

In  order  to  make  the  computational  grid  for  the  FD  method  finite,  it  is  necessary 
to  impose  some  boundary  conditions  on  the  elastic  field  at  finite  values  of  x  and  2^. 
The  true  boundary  conditions  are  that  the  scattered  portions  of  the  field  should  be 
only  outgoing  and  it  is  possible  to  simulate  these  boundary  conditions  with  varying 
degrees  of  accuracy  [13].  The  approach  of  the  FD  code  of  this  report  is  to  set  the 
elastic  variables  to  zero  at  the  edges  of  the  grid.  This  causes  a  spurious  reflection 
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of  energy  incident  on  the  sides  of  the  grid.  In  order  to  eliminate  this  artificially 
reflected  energy,  an  attenuating  layer  is  introduced  to  the  grid,  running  parallel  to 
the  sides  of  the  grid.  Two  MATLAB  modules  have  been  written  for  the  creation 
of  this  attenuating  layer.  The  first  module  foUows  the  suggestion  of  [5]  and  allows 
the  user  to  define  a  very  low  Q  value  at  the  edge  of  the  grid  (e.g.  Q  ~  2),  A 
smooth  transition  from  the  true  Q-value  in  the  interior  of  the  grid  to  the  low  value 
at  the  edge  is  used  (typically,  a  transition  over  20-40  grid  points  is  used).  Thus  we 
are  applying  the  visco-elastic  capabilities  of  the  code  to  the  problems  of  artificial 
reflections  from  the  grid  boundaries.  As  was  explained  in  the  previous  subsection 
(see  Fig. 2)  there  is  significant  dispersion  for  low  values  of  Q.  In  order  to  eliminate 
reflections  caused  by  the  change  of  phase  velocity  due  to  changing  Q,  the  compres- 
sional  and  shear  velocities  defined  on  the  grid  are  adjusted  using  Eq.(lO)  so  that  the 
phase  velocity  at  the  dominant  frequency  of  the  problem  is  a  constant  with  respect 
to  Q.  Automatic  velocity-tuning  is  a  feature  of  the  MATLAB  module.  A  second 


Figure  5:  Attenuating  factor  from  Eq.(16)  for  a  transition  layer  of  fO  grid  points 
using  p=2  (solid  line)  and  p=0.4  (dashed  line). 


attenuation  mechanism  is  implemented  by  multiplying  the  wavefield  within  the  at¬ 
tenuating  layer  solution  by  a  constant  (  at  each  time  step.  The  constant  (  is  equal 
to  one  at  the  start  of  the  attenuating  layer  (i.e.,  within  the  interior  where  the  layer 
starts  and  is  tapered  to  a  smaller  value  (e.g.  =  0.97)  at  the  edge  of  the  grid. 

The  MATLAB  module  which  was  written  for  defining  these  layers,  uses  a  tapering 
function  of  the  form 


C(i)  =  (1  -  Uin) 


1  +  cos(S 


+  Cv 


i  = 


(16) 
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This  type  of  approach,  in  conjunction  with  radiating  boundary  conditions  was  used 
by  Fricke[l],  In  practice,  when  using  Eq.(16),  we  used  p=:2  or  p=0.4  and  N=  20 
or  40  points.  In  Fig. 5  we  show  the  attenuating  function  of  Eq.(16)  (with  Cmm  — 
0.97)  for  N=40  and  p=2  (solid  line)  and  p=0.4  (dashed  line).  The  two  curves  are 
quite  different  in  character.  The  p  =  2  curve  has  its  largest  slope  in  the  middle  of 
the  transition  zone  whereas  the  p  =  0.4  curve  is  gradually  decaying  for  much  of  the 
curve  and  then  has  a  large  (infinite  at  the  end  point)  slope  at  the  end.  We  found 
that  for  modelling  a  point  source  the  p  —  0.4  curve  worked  well;  however,  in  general, 
one  must  often  experiment  with  the  absorbing  boundary  in  order  to  ensure  that  the 
FD  solutions  are  sufficiently  free  of  spurious  boundary  reflections. 

For  the  FD  code  of  this  report,  one  can  use  either  of  the  attenuating  mechanisms 
described  above  or  both.  For  the  computations  of  this  report  we  used  only  the 
second  technique. 


2.4  Specification  of  incident  field 

In  our  implementation  of  the  FD  code  we  use  point  excitations  at  grid  points  to 
generate  incident  fields.  We  consider  below  some  of  the  details  of  modelling  a  point 
source  which  lies  within  the  boundaries  of  the  numerical  grid.  Second,  we  discuss 
the  generation  of  an  incident  field  from  a  distant  source  by  using  an  array  of  point 
sources. 


Interior  point  source 


The  basic  equations  of  elastic  propagation  (ignoring  any  additional  visco-elastic 
terms)  can  be  written  as 


Vx,t 


^xx,t 


^zz,t 


^xZfi 


I  f  ^^XX  J  ^^xz'^. 

p\  dx  dz  ) 

1  /  da^z  dcTzz\ 
p\  dx  dz  ) 


(A  +  2p,)-^  +  A 


dvz 

dz 

dvx 

dx 


(17) 

(18) 

(19) 

(20) 
(21) 


Let  us  consider,  in  an  acoustic  medium,  a  point  source  located  at  a;  =  (discrete 
grid  point  i^)  and  z  =  Zg  (discrete  grid  point  jg)  which  emits  a  signal  S{t).  Since  the 
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medium  is  acoustic,  the  terms  /i  and  can  be  set  to  zero  in  the  above  equations. 
Also,  (Jxx  =  ^zz  =  Equation(20)  can  now  be  written  as 

(iv  3v 

at  =  +  S{t)6{x  -  Xs)S{z  -  2s).  (22) 

Differentiating  Eq.(22)  with  respect  to  time  and  utilizing  Eqs.(17)  and  (18)  we  obtain 


di'^  p  \dx^  dz^ ) 


+ 


dS{t) 

dt 


S{x  -  Xs)S{z  -  Zs) 


(23) 


or 

~  =  c^V'^a  +  St{t)S{x  -  Xs)S{z  -  Zs).  (24) 

In  the  discretized  version  of  Eq.(20)  the  source  term  is  implemented  using  the  equa¬ 
tion, 


^Xxi'^Sf  js) 


^ZziSs")  js) 


^xx{,is  js)  H" 
^zzi^is'i  js)  ■!" 


S{t){At) 

AxAz 

(25) 

S{t){At) 

(26) 

The  normalization  with  respect  to  At  is  due  to  the  time- discretization  of  Eqs.(19) 
and  (20);  the  spatial  normalization  is  chosen  so  that  the  source  appears  as  a  spatial 
delta- function  with  respect  to  discrete  integration  [14].  From  Eq.(24)  it  can  be  seen, 
that  in  order  to  compare  the  ED  solutions  with  analytical  solutions  of  the  acoustic 
wave  equation,  it  is  necessary  to  use  the  time-derivative  of  the  source  function  S{t) 
as  the  source  function  in  the  analytical  solution. 


An  incident  field  from  a  distant  source 

Because  it  is  necessary  to  use  10-20  grid  points  per  wavelength  in  the  ED  mod¬ 
elling,  the  ED  grid  corresponds  to  relatively  small  physical  dimensions  for  higher 
frequencies.  However,  for  many  problems  of  interest  the  source  may  be  so  far  away 
from  the  scattering  object  or  surface,  that  it  is  not  feasible  to  include  the  source 
point  within  the  computational  grid.  If  the  waveguide  between  the  source  and  the 
scattering  region  is  sufficiently  simple,  it  may  be  possible  to  analytically  compute 
the  incident  field  at  the  edges  of  a  numerical  grid.  From  the  integral  relations  of 
the  Gauss’s  divergence  theorem  [15]  we  can  use  these  boundary  values  to  excite  the 
ED  code  and  propagate  the  incident  field  within  the  grid.  This  type  of  approach 
has  been  used  previously  by  authors  in  electromagnetic  modelling  (see  for  example, 
[16]  where  an  approach  very  similar  to  that  considered  in  this  report  is  used).  We 
now  present  some  details  of  our  particular  implementation.  We  consider  the  case 
of  an  acoustic  waveguide  between  the  source  and  scattering  region.  From  Gauss’s 
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Divergence  Theorem,  we  can  express  the  acoustic  pressure  field  at  a  single  frequency 
ij  as 

p(^)  =  -  Gn’p)di  (27) 

where  C  is  a  bounding  rectangle  located  within  the  computational  grid,  G  is  a  waveg¬ 
uide  Green’s  function  satisfying  the  Helmholtz  equation 

=  6{x  -  x')6{z  -  z'),  (28) 

c^{z) 

and  n'  denotes  the  normal  derivative  with  respect  to  the  coordinates  of  integration 
along  the  boundary.  We  can  interpret  the  integral  of  Eq.(27)  as  the  pressure  field 
resulting  from  a  distribution  of  monopole  and  dipole  sources  (Huygen’s  sources  [17]). 
If  the  source  is  distant  from  the  scattering  region,  then  the  incident  field  will  look 
like  an  incident  plane- wave  (or  as  we  wiU  discuss  a  generalized  plane- wave  for  a  two 
halfspace  medium)  in  the  water  column  and  this  plane- wave  wiU  have  an  associated 
angle  of  incidence.  We  define  the  numerical  grid  so  that  one  of  its  edges  is  normal 
to  the  angle  of  incidence  of  the  plane- wave  (Fig.  6)  and  is  internal  to  the  absorbing 
boundary  layer  of  the  grid.  We  consider  Eq.(27)  only  along  this  line  and  ignore 
the  other  boundaries’  contributions  to  the  incident  field.  There  are  different  ways 
of  implementing  Eq.(27).  We  can  consider  a  line  of  appropriately  weighted  discrete 
monopole  point  sources  along  the  boundary.  The  ED  solution  for  this  array  of 
sources  effectively  performs  the  integration  with  respect  to  the  Green’s  function  ( 
a  convolution  in  the  time  domain).  Because  of  the  symmetrical  field  produced  by 
these  monopole  sources,  Gx'  is  zero  along  the  line  and  we  only  have  the  first  term 
in  the  integral  of  Eq.(27).  Alternatively,  we  could  simulate  a  series  of  dipole  sources 
by  placing  a  sequence  of  (-h)  and  (-)  monopoles  one  grid  point  on  either  side  of  the 
bounding  edge.  In  this  case  only  the  second  term  of  the  integral  is  required.  A 
simpler  method  to  implement  a  dipole  is  to  add  a  spatial  delta  function  (weighted 
by  the  incident  pressure  field  normalized  by  density)  to  the  equation,  Eq.(17),  for 
Vx’  For  aU  these  boundary  implementations,  the  sources  produce  a  field  which 
propagates  outwards  in  both  directions.  Thus  not  only  the  incident  field  is  produced, 
but  also,  a  field  incident  on  the  grid  boundary.  If,  however,  the  attenuating  layer 
is  working  weU,  this  field  will  not  reenter  the  grid.  If  Eq.(27)  is  implemented  using 
both  the  monopoles  and  dipoles  then  we  can  produce  only  the  incident  field.  The 
approach  implemented  in  the  FD  code  was  to  excite  the  equation  for  Vx{z)  along 
the  line  x  —  Xinc  with  P{xiyic^  p{^)^ 

At  this  point  our  discussion  has  been  general  in  terms  of  form  of  the  incident  field; 
however,  we  have  numerically  implemented  the  case  where  the  incident  pressure  field 
is  a  plane  wave  with  an  angle  of  incidence  9i  in  a  two  halfspace  environment.  Then 
the  incident  field  consists  of  the  direct  wavefield 

P^{x^z,t)  =  Reall  f  S{u)exp  iuj{t  -  (—cos(0i)  — —sin(6i))  doA^  (29) 
Uo  L  Cl  Cl  J  J 
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the  reflected  wavefield 


t)  =  Real 


I  /  S{u})R{9i,uj)exp  iuj{t  —  {—cos{6i) -\ - sin{6i))  duf,  (30) 

Uo  [Cl  Cl  J  J 


and  a  transmitted  field, 


P^{x,z,t)  =  Really  j  S{uj)T{6i,uj)exp 


iult  —  i — cos{0t) - siniOt))  du)\,  (31) 

C2  C2  J  J 


where 

P2C2sin{ei)  -Pii/cj-  clcos^(ei) 

Ryui^uj)  —  I - 5  (32) 

P2C2sin{6i)  -\-  piijc\-  c\cos^{6i) 
p2C2sin{di)  +  P\yc\-  clcos^{9i) 

and  S{lo)  is  the  Fourier  Tranform  of  the  source  function.  The  factor  yjc\  —  c\cos^{6i) 
is  imaginary  for  Oi  less  than  the  critical  angle,  in  which  case  the  reflection  and 
transmission  coefficients  are  complex.  This  means  that  in  the  time  domain  the 
reflected  and  transmitted  pulses  are  combinations  of  the  source  pulse  shape  and  its 
Hilbert  Transform  [18]. 


Pi^i 


Figure  6:  Schematic  diagram  of  orientation  of  Finite  Difference  grid  with  respect  to 
a  distant  source 
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3 

Numerical  Examples 


In  this  section  we  first  show  the  numerical  results  obtained  from  the  finite  differ¬ 
ence  modelling  of  propagation  and  scattering  problems,  for  which  we  can  also  com¬ 
pute  analytical  solutions.  In  particular,  we  consider  the  propagation  of  an  acoustic 
wavefield  from  a  point  source  in  non-attenuating  and  attenuating  media.  We  then 
consider  the  propagation  of  a  shear  plane- wave  in  an  attenuating  medium.  Finally, 
for  the  benchmark  cases,  we  consider  the  scattering  of  a  compressional  plane-wave 
by  solid,  thick-shelled,  and  thin-sheUed  aluminum  cylinders.  The  shelled  cylinders 
are  filled  with  water.  The  analytical  solutions  in  this  case  are  computed  by  solving 
in  the  frequency  domain  a  sequence  of  elastic  cylinder  scattering  problems  and  then 
constructing  the  time-domain  pulse  by  Fourier  synthesis. 

After  these  benchmark  cases,  we  consider  a  buried  solid  aluminum  cylinder  and  con¬ 
sider  a  generalized  (i.e.,  consisting  of  direct,  reflected,  and  transmitted  components) 
plane- wave  incident  upon  this  object.  We  consider  the  cases  of  the  direct  wave  in 
the  water  column  having  pre-critical  and  post-critical  angles  of  incidence.  For  the 
case  of  an  evanescent  transmitted  wave,  we  examine  the  effect  of  a  rough  interface. 
Finally  we  use  the  FD  method  to  compute  the  backscattered  field  from  a  cylinder 
at  an  array  of  receivers  in  the  water  column  and  then  time-reverse  this  field  and  use 
it  as  the  “incident”  excitation  field.  The  resulting  backpropagated  field  in  the  water 
column  focuses  at  the  sources  of  scattering. 

The  numerical  code  was  run  in  Fortran  on  a  DEC-3000  Alpha  workstation.  A 
computation  with  a  440  X  440  grid  and  12000  time  steps  required  128  minutes  of 
CPU  time.  In  the  following  examples,  the  absolute  levels  of  the  computed  pressure 
fields  have  often  been  scaled  for  plotting  purposes. 


3.1  Computations  in  a  homogeneous  medium 


First  we  consider  a  point  source  in  a  homogeneous  acoustic  space.  We  consider  the 
source  function  in  Eq.(22)  to  be  the  time  derivative  of  a  Gaussian  pulse,  [1], 

5(0  = (34) 
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where 


a  = 


(35) 

(36) 


and  cjo  is  the  central  angular  frequency  which  we  take  to  be  1500  Hz.  We  take  the 
sound  speed  to  be  cq  —  1500  m/s.  Thus  the  reference  wavelength  for  this  example 
is  Im.  We  do  our  computations  with  Ax  =  0.05  m  or  A/20.  From  our  previous 
discussion  of  stability  limits,  we  must  use  a  time  step  such  that 


At  <  .606 


Ax 

1^ 


2  X  10”^ 


s. 


(37) 


We  use  =  10”^5ec  which  is  approximately  half  the  stability  bmit.  We  consider  a 
numerical  grid  which  is  300  X  300  in  dimension.  We  use  40  points  for  the  absorbing 
boundary  layer.  We  construct  this  absorbing  layer  using  the  second  absorbing  mech¬ 
anism  described  in  subsection  2.4.  We  found  that  using  a  fractional  power  p  =  0.4 
in  Eq.(16)  with  (  =  0.95  worked  well.  In  Fig. 7  below,  we  show  the  FD  signal  as  a 
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Figure  7:  Finite  difference  (solid  line)  and  analytical  (dashed  line)  pulse  shapes  for 
receiver  at  8.5  wavelengths  from  the  source 

solid  line  along  with  the  analytical  solution  of  Eq.(24)  (dashed  line)  for  a  source  at 
{i,j)  =  (70, 150)  and  a  receiver  at  {i^j)  =  (240, 150)  (where  (i^j)  denote  the  discrete 
indices  for  the  (x,  z)  coordinates).  The  agreement  is  excellent  between  the  computed 
and  analytical  solutions  with  only  small  artifacts  due  to  the  boundaries  in  the  tail 
of  the  signal. 
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We  now  consider  the  same  geometry  as  above,  but  for  a  medium  with  Q  =  20.  In 
Figure  8a  we  show  the  comparison  between  the  FD  computed  pulse  and  an  “an¬ 
alytically”  computed  pulse.  The  analytical  solution  used  Fourier  synthesis  and  a 
Q/ dispersion  relation  of  the  form  of  Eqs.(14)  and  (15)  with  the  parameters  which 
were  used  for  the  computation  of  the  curve  of  Fig.  3.  The  agreement  is  excellent. 
In  Fig.  8b  we  show  a  comparison  between  the  FD  pulse  and  the  analytical  solution 
if  no  frequency  dispersion  is  used;  i.e.,  we  fixed  the  phase  velocity  at  c  =  1500  m/s 
for  aU  frequencies  and  simply  added  an  imaginary  part  to  the  sound  speed  to  pro¬ 
duce  the  required  attenuation.  This  procedure  is  the  usual  approach  in  frequency 
domain  modelling.  The  agreement  between  the  two  pulses  is  no  longer  so  good.  The 
amplitudes  of  the  two  pulses  are  similar  but  the  shape  of  the  pulses  are  noticeably 
different.  We  have  used  a  large  attenuating  factor  in  this  example;  the  differences 
between  the  curves  in  Fig.  8b  would  be  less  for  smaller  attenuation.  Finally  in  Fig.  8c 
we  show  the  analytical  pulse  computed  for  a  Q=1000  (no  attenuation).  In  this  case 
the  pulse  is  too  large  in  amplitude,  illustrating  the  importance  of  attenuation  in  this 
example. 


Time(ms) 


Figure  8:  Finite  difference  (solid  line)  and  analytical  pulse  shapes  (dashed  line)  for 
Q  =20  for  (a)  analytical  solution  computed  accounting  for  frequency  dispersion  (b) 
attenuation  model  with  a  fixed  phase  velocity  oflSOOm/s  (c)  using  Q=1000 


We  now  consider  an  example  with  a  non-zero  shear  speed.  In  particular  we  excite 
the  equation  for  Vz  with  an  array  of  point  sources  along  the  line  x  =  70.  This 
arrangement  should  excite  a  shear  plane-wave.  However,  due  to  the  fact  that  there 
is  a  finite,  discrete  aperture  of  sources  the  synthesis  of  the  plane-wave  is  not  exact 
and  other  energy,  some  of  it  compressional  in  nature,  is  excited.  We  use  an  absorbing 
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layer  of  40  points  in  this  case,  with  (  =  0.96  and  p  =  2.  To  reduce  the  amount  of 
spurious  energy  produced  by  the  edge  effects  of  a  finite  aperture  of  sources,  we  weight 
the  source  terms  by  a  cosine  taper  away  from  the  centre  of  the  grid.  The  source 
time  history  is  the  same  as  for  the  previous  example.  However,  we  now  take  the 
compressional  sound  speed  to  be  Cp  =  3000m/5  and  the  shear  speed  —  1500m/ s. 
Because  of  the  high  compressional  sound  speed  we  use  At  =  10”^ /2.  We  define  the 
Q  factor  for  the  compressional  waves  to  be  Qp  =  1000  throughout  the  grid  and  the  Q 
factor  for  the  shear  waves  =  20.  If  the  FD  code  is  correctly  modelling  the  elastic 
and  attenuation  properties  of  the  medium,  the  wave  propagation  for  this  example 
should  be  weU-modeUed  as  a  shear  plane  wave  attenuated  by  a  Q  value  of  20.  In 
Fig. 9  we  show  the  pulse  (we  plot  v^)  at  the  centre  of  the  grid  at  a:  ~  240  as  computed 
by  the  FD  code  (solid  line)  and  analytically  (dashed)  for  an  attenuated  shear  plane 
wave.  The  agreement  between  the  two  curves  is  good.  There  is  a  small  amount  of 
spurious  energy,  some  of  which  is  compressional  (note  that  any  compressional  energy 
is  much  less  attenuated  than  the  true  shear  energy  as  we  have  used  Qp  =  1000  for 
this  example).  It  is  interesting  to  note  that  for  this  example,  the  Green’s  function 
which  we  use  in  the  computation  of  the  analytic  solution  is  of  the  form 

exp{iuj/c\x  —  a^sl) 
c 

whereas  in  the  point  source  example  the  Green’s  function  was  of  the  form 

jH^{u/cy/{x  -  XsY  +  {z-  Zs)^)  (39) 


(38) 


3.2  Computations  for  an  elastic  cylinder  surrounded  by  an  homogeneous  fluid 

We  now  consider  an  aluminum  cylinder  of  radius  1  m  located  in  an  homogeneous 
compressional  space  with  Cp  =  1500m /s  and  p  ~  lOOOkg/m^.  For  the  aluminum 
cylinder  we  use  the  parameters  Cp  —  6380m/5,  Cs  =  31367n/5  and  p  =  2172%/m^. 
Because  of  the  large  velocities  within  the  cylinder  it  is  necessary  to  take  small  time 
steps.  We  used  Ax  ^  0.025m  (or  40  points  per  wavelength)  and  At  —  0.01/8m5  in 
our  computations  with  a  400  X  400  grid.  The  cylinder  is  located  at  the  centre  of 
the  numerical  grid.  A  normal-incidence  plane-wave  was  excited  by  applying  point 
source  functions  to  and  a^z  along  the  line  x  =  1.5m  with  the  time-history  of  a 
time-differentiated  gaussian.  In  order  to  reduce  edge  effects  we  tapered  the  incident 
plane-wave  by  a  cosine- weighting  for  the  first  and  last  100  points  of  the  grid.  Figure 
10  shows  a  comparison  of  the  FD  time  series  with  the  time  series  computed  using 
Fourier  synthesis  and  the  harmonic  Fourier-Bessel  series  for  the  scattered  field  from 
an  infinite  elastic  cylinder.  The  series  were  computed  for  receivers  3.5  m  from  the 
centre  of  the  cylinder  at  angles  of  0"^,  90^,  270^  and  180®  with  respect  to  the  cylinder 
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Figure  9;  Finite  difference  (solid  line)  and  analytical  pulse  shapes  (dashed  line)  for 
a  shear  plane-wave  with  Q  =  20 


centre.  From  symmetry,  the  time  series  for  9  =  90°  and  270°  should  be  identical. 
However,  this  is  not  true  in  practice  due  to  the  staggered  grid  formulation  of  the  RBS 
FD  code.  Although  the  FD  solution  is  initially  symmetrical,  the  symmetry  of  the 
field  is  not  sustained  over  time.  The  FD  and  analytical  amplitudes  of  the  incident 
wave  do  not  agree  well  for  the  receivers  at  90°  and  270°  because  the  tapering  of 
the  plane-wave  is  significant  at  these  locations.  The  scattered  fields,  however,  are 
caused  by  the  scattering  of  the  incident  field  at  the  cylinder  and  these  amplitudes 
are  correctly  modelled. 

We  now  repeat  the  above  numerical  experiment  with  a  shelled  cylinder.  We  consider 
an  aluminum  shell  0.5  m  thick  (or  20  grid  points)  with  the  interior  of  the  cylinder 
filled  with  water.  The  comparison  with  the  analytic  solution  is  shown  in  Fig.  11. 
The  agreement  between  the  numerical  and  analytic  solutions  is  good. 

By  varying  the  shell  thickness,  we  found  that  good  agreement  between  computed 
and  analytical  results  were  obtained  down  to  and  including  shell  thicknesses  of  3 
grid  points  corresponding  to  a  relative  shell  thickness  of  7.5%;  by  decreasing  the 
spatial  discretization  size  we  should  be  able  to  decrease  this  ratio.  The  comparison 
between  the  FD  curves  and  the  analytical  curves  are  shown  in  Fig.  12.  In  Fig.  13 
we  show  a  snapshot  of  the  pressure  field  for  the  thick  shelled  cylinder  4  ms  into  the 
computation.  The  plane-wave  can  be  seen  passing  by  the  exterior  of  the  cylinder. 
A  faster  wavefront  has  already  passed  through  the  cylinder  and  the  backscattered 
wavefront  is  visible  at  the  front  of  the  cylinder. 
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Figure  10:  Finite  difference  (solid  line)  and  analytical  time  series  (dashed  line)  for 
receivers  at  (a)  0°  (backscatter)  (b)  90°  (c)  270°  and  (d)  180°  (forward  direction) 
for  an  infinite  elastic  (aluminum)  cylinder 


Figure  11:  Finite  difference  (solid  line)  and  analytical  time  series  (dashed  line)  for 
receivers  at  (a)  0°  (backscatter)  (b)  90°  (c)  270°  and  (d)  180°  (forward  direction) 
for  thick-shelled  (20  grid  points)  cylinder  filled  with  water 
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Figure  12:  Finite  difference  (solid  line)  and  analytical  time  series  (dashed  line)  for 
receivers  at  (a)  0°  (backscatter)  (b)  90^  (c)  270°  and  (d)  180°  (forward  direction) 
for  thin-shelled  (3  grid  points)  cylinder  filled  with  water 


3.3  Scattering  of  a  generalized  plane-wave  by  a  buried  cylinder 

We  can  simulate  an  incident  acoustic  or  elastic  wavefield  by  exciting  monopole  and 
dipole  sources  along  the  boundaries  of  the  grid.  In  these  examples,  we  only  use 
dipole  excitation.  We  excite  the  particle  velocity  by  p^’^^ j p{z).  This  means  that 
the  correct  field  propagates  into  our  grid,  but  also  a  mirror-image  field  propagates 
to  the  left  where  it  is  attenuated  by  the  absorbing  layer  of  our  grid.  By  exciting  the 
normal  stress  by  it  should  be  possible  to  eliminate  the  left-going  field. 

In  the  examples  which  foUow,  we  wiU  consider  a  fluid  half-space  with  sound  speed 
1500  m/s  and  p  =  lOOO/c^f/m^  overlying  another  fluid  space  with  sound  speed  1700 
m/s  and  p  =  1500A:5r/m^.  The  critical  angle  for  these  parameters  is  9^  —  28.07° 
(measured  from  the  horizontal).  We  first  consider  the  case  of  the  angle  of  incidence 
9  —  30°.  We  consider  our  grid  rotated  so  that  the  boundary  of  excitation  is  parallel 
to  the  incident  direct  wave  and  as  a  result  the  fluid/fluid  interface  has  a  slope  of  30°. 
The  advantage  of  this  approach  is  that  one  obtains  a  maximum  aperture  of  sources, 
using  a  single  boundary,  to  simulate  an  incident  field.  The  disadvantages  are  (1) 
that  the  interface  is  sloping  upwards  to  the  right  of  the  incident  excitation  line  and 
it  is  difficult  to  know  how  to  model  the  bathymetry  to  the  left  of  the  incident  line 
in  order  to  avoid  diffraction  effects  at  the  intersection  of  the  interface  boundary  and 
the  vertical  array  of  excitation  sources;  (2)  because  the  interface  is  sloping  we  are, 
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Figure  13:  Snapshot  of  pressure  field  at  4  ms  for  thick-shelled  aluminum  cylinder 
filled  with  water 
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in  fact,  approximating  it  with  a  staircase  in  the  FD  approximation  and  these  small 
stairsteps  cause  small  diffractions.  However,  we  find  that  despite  these  problems 
the  incoming  incident  field  is  well-behaved.  In  Fig.  14  we  show  the  incident  field 
(i.e.,  'p{z^t)l p[z))  which  is  used  to  excite  The  vertical  axis  is  depth  and  the 
horizontal  axis  is  time.  The  direct  pulse  appears  vertically  oriented.  As  can  be  seen 
there  is  significant  energy  transmitted  into  the  bottom.  For  the  FD  computations 
of  this  set  of  examples,  we  used  a  440  X  440  grid  with  60  points  in  the  absorbing 
layer  (Cmm  =  0.96,  p  =  2).  The  equations  for  were  excited  at  the  line  Xinc  —  63. 
In  Fig.  15  we  show  a  time  snapshot  of  the  pressure  field  after  the  plane  wave  in  the 
bottom  has  interacted  with  the  buried  aluminum  cylinder.  There  is  a  substantial 
reflected  field  from  the  cylinder.  Only  the  field  in  the  square  region  interior  to  the 
absorbing  boundary  layers  and  the  line  of  dipoles  is  shown. 
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Figure  14:  Incident  plane-wave  field  (absolute  value)  used  to  excite  vertical  array  of 
dipole  sources  for  6  =  30°. 


Next,  we  consider  the  case  for  the  angle  of  incidence  of  the  direct  plane- wave  in 
the  water  column  being  post-critical,  9  —  20°.  The  incident  field  used  to  excite  the 
equation  for  Vx  is  shown  in  Fig.  16. 

In  this  case  the  transmitted  field  is  evanescent  and  the  reflected  pulse  is  a  temporal 
combination  of  the  direct  pulse  and  its  Hilbert  Transform.  There  is  a  region  of  high 
intensity  near  the  interface  where  the  incident  and  reflected  waves  interfere  with 
each  other.  In  Figs.  17,  18,  and  19  we  show  a  single  time  snapshot  for  the  wavefields 
generated  in  the  case  of  a  buried  solid  aluminum  cylinder,  for  a  rough  interface  with 
no  cylinder,  and  the  same  rough  interface  with  a  buried  cylinder.  In  Fig.  17  we  see 
that  the  incident  field  of  Fig.  16  has  propagated  into  the  grid.  The  exponential  tail 
of  the  transmitted  wave  just  touches  the  cylinder.  There  is  a  small  amount  of  energy 
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Figure  15:  Time  snapshot  of  pressure  field  (absolute  value)  resulting  from  interaction 
of  transmitted  wavefield  with  cylinder 
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Figure  16:  Incident  plane-wave  field  (absolute  value)  used  to  excite  vertical  array  of 
dipole  sources  for  6  =  20° . 
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Figure  17:  Propagating  plane-wave  field  (absolute  value)  with  buried  aluminum  cylin¬ 
der 
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Figure  18:  Propagating  plane-wave  field  (absolute  value)  with  rough  interface 
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Figure  20:  Time  series  at  a  receiver  located  3.25  m  above  interface  at  incident  array 
when  there  is  no  cylinder  in  the  bottom  (solid  line)  and  when  there  is  a  cylinder  in 
the  bottom  (dashed  line) 
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following  the  incident  wavefield  which  was  probably  incorrectly  created  during  the 
generation  of  the  incident  wavefield  from  the  dipoles.  This  can  sometimes  be  an 
annoyance  when  generating  incident  post-critical  wavefields;  it  is  not  surprising  that 
the  use  of  a  discontinuous  (because  of  the  p{z)  normalization),  analytic  wavefield 
to  excite  a  discrete  numerical  grid  may  cause  the  excitation  of  a  small  amount  of 
spurious  energy  near  the  interface.  It  may  be  possible  to  reduce  this  problem  by 
smoothing  the  incident  wavefield  and/or  the  medium,  but  we  did  not  do  this  for 
these  results. 


The  “rough”  interface  of  Figs.  18  and  19  is  generated  from  the  equation 

j  =  170  +  (i  -  63)  tan(20°)  +  10  (40) 


Here  j  is  the  vertical  index  and  i  the  horizontal  index.  There  are  40  points  per  wave¬ 
length  used  in  this  example  so  that  the  peak-to-peak  roughness  is  1/2  wavelength 
and  the  wavelength  of  the  roughness  is  one  wavelength.  The  rough  interface  gen¬ 
erates  backscattered  energy  and  also  transmits  energy  into  the  bottom  (Fig.  18). 
This  energy  can  significantly  interact  with  the  cylinder  (Fig.  19)  -  however,  the 
differences  between  Figs.  18  and  19  in  the  water  column  are  very  slight.  This  is 
emphasized  in  Fig.  20  where  we  show  the  received  timeseries  for  a  receiver  3.25  m 
above  the  interface  at  the  horizontal  location  of  the  array  of  Huygen’s  sources  for  the 
case  of  no  cylinder  in  the  bottom  (solid  line)  and  a  cylinder  in  the  bottom  (dashed 
line).  The  two  curves  are  almost  indistinguishable.  There  is,  however,  a  small  but 
noticeable  difference  between  the  two  curves  in  the  6-8  ms  interval.  Another  inter¬ 
esting  feature  of  these  two  curves  is  the  the  appearance  of  a  definite  frequency  of 
backscatter  for  this  case.  This  is  expected  from  a  perturbation  analysis  of  the  rough 
interface.  Following  the  work  of,  for  example,  [19]  we  expect  for  our  given  roughness 
wavelength  that  the  maximum  backscattered  field  should  occur  approximately  when 


2 


u 


=  27rA(perturbation) 


(41) 


In  our  case  the  perturbation  has  a  wavelength  of  1  m  so  that  we  find  an  optimal 
frequency  of  approximately  750  Hz  or  a  period  of  oscillation  of  1.3  ms  which  agrees 
weU  with  the  result  of  Fig.  20. 


3.4  An  example  of  time  reversed  propagation 

As  shown  above  it  is  straightforward  to  introduce  an  incident  field  into  the  grid 
by  using  an  array  of  monopoles  and/or  dipoles.  Another  interesting  application  of 
this  concept,  is  to  time- reverse  [20]  the  scattered  field  at  the  receiving  array  and 
propagate  it  back  into  the  medium.  We  would  expect  the  field  which  has  been 
scattered  from  compact  objects  to  focus  back  at  those  same  objects. 
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We  consider  the  same  cylinder  and  medium  as  the  previous  example  but  now  the 
cylinder  is  only  partially  buried.  A  time  snapshot  of  the  pressure  field  is  shown  in 
Fig.  21  just  after  the  incident  field  has  encountered  the  cylinder.  In  Fig.  22  we 


ABOVE  1.05 
0.85-1.05 
0.65  -  0.85 
0.45  -  0.65 
0.25  -  0.45 
0.05  -  0.25 
BELOW  0.05 


Figure  21:  Plane-wave  field  incident  on  partially  buried  aluminum  cylinder 

show  the  time-reversed  pressure  field  recorded  by  210  receivers.  The  first  of  these 
receivers  is  located  in  the  water  just  above  the  interface  and  the  others  are  located 
sequentially  above  it.  Starting  from  the  right  hand  side  of  the  plot,  the  incident  field 
is  evident,  then  backscatter  from  the  rough  interface,  and  then  from  7.0  to  5.5  ms 
an  event  associated  with  backscatter  from  the  cylinder,  and  then  finally  backscatter 
from  the  rough  interface  again.  The  field  of  Fig.  22  is  then  used  to  excite  the  array  of 
dipoles  (just  for  the  positions  of  the  210  receivers).  We  stop  the  time-reversed  field 
just  prior  to  the  time  that  the  original  incident  field  is  present  in  the  time  series.  In 
this  numerical  experiment  of  time-reversed  propagation,  we  suppose  that  we  do  not 
know  what  created  the  scattered  field.  Hence  although  the  cylinder  was  used  in  the 
FD  modelling  to  produce  the  scattered  field  of  Fig.  22,  it  is  not  included  in  the  model 
for  the  backpropagation  of  this  energy.  Instead,  only  a  two  halfspace  medium  is  used 
in  the  modelling  of  the  backpropagated  scattered  field.  In  Figs.  23  and  24  we  show 
2  time  snapshots  of  the  resulting  field  as  it  propagates  back  into  the  grid.  Figure  23 
shows  the  field  scattered  from  the  cylinder  just  starting  to  focus  and  Fig.  24  shows 
the  field  close  to  the  time  of  maximum  focussing;  the  focus  is  along  the  side  where 
the  original  specular  reflection  occurred.  For  this  example,  our  recording  aperture 
did  not  sufficiently  surround  the  scatterer  and  the  frequency  of  the  source  was  not 
high  enough  to  obtain  a  good  definition  of  the  scattering  object;  nevertheless,  the 
backpropagated  field  did  focus  at  the  correct  location,  corresponding  to  the  area  of 
specular  reflection  from  the  object. 
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Figure  22:  The  time  history  of  the  time-reversed  pressure  field  (absolute  value)  as 
recorded  along  a  vertical  array  of  receivers  located  at  the  line  of  incidence 
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Figure  23:  Back-propagated  field  (absolute  value)  as  scattered  field  is  just  starting  to 
focus;  the  cylinder  is  not  used  in  the  modelling 
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Figure  24:  Back-propagated  field  (absolute  value)  near  the  time  of  maximum  focus¬ 
ing;  the  cylinder  is  not  used  in  the  modelling 
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Summary 


We  have  described  some  of  the  theory  underlying  a  visco-elastic  Finite  Diiference 
code  (in  particular,  the  RBS  FD  code).  We  compared  the  FD  computations  with  a 
suite  of  analytical  cases  for  both  homogeneous  propagation  and  scattering  and  found 
that  the  FD  results  were  in  good  agreement  with  the  analytical  results.  The  attenu¬ 
ation  properties  of  the  FD  code  are  in  good  agreement  with  analytical  computations 
if  one  uses  a  realistic  frequency  dipersion  relation  in  the  analytical  modelling.  We  did 
computations  of  scattering  from  aluminum  cylinders,  which  have  a  very  high  veloc¬ 
ity/density  contrast  with  the  surrounding  fluid  and  obtained  good  agreement  with 
analytical  computations  for  the  4  quadrants  of  scattering,  although  the  results  were 
best  for  the  backscatter  direction.  We  consider  aluminum  shells  of  varying  thickness 
and  found  that  we  obtained  good  results  for  shells  only  3  grid  points  thick. 

We  then  showed  how  we  could  introduce  an  incident  field  into  the  grid  by  exciting 
an  array  of  dipoles.  We  used  this  concept  to  introduce  a  generalized  plane-wave, 
corresponding  to  a  two  half-space  medium  into  the  grid.  In  particular  we  considered 
an  incident  wave  for  which  the  wavefield  in  the  bottom  half-space  was  evanescent 
and  considered  the  scattering  by  a  buried  aluminum  cylinder.  We  then  repeated  the 
computations  for  a  rough  interface.  For  this  particular  example,  the  rough  interface 
introduced  more  energy  into  the  bottom  but  the  signal  received  in  the  water  column 
was  dominated  by  the  backscatter  from  the  rough  interface.  Reflected  energy  from 
the  cylinder  was  just  observable  in  the  backscattered  signal.  Finally,  we  showed 
that  with  the  incident-field  formulation  of  this  code,  one  can  easily  back-propagate 
scattered  fields  and  focus  this  energy  into  areas  of  high  reflectivity. 
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